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Abstract 

We study the evolution of fronts in a bistable reaction-diffusion system when the 
nonlinear reaction term is spatially non- homogeneous. This equation has been used 
to model wave propagation in various biological systems. Extending previous works 
on homogeneous reaction terms, we derive asymptotically an equation governing the 
front motion, which is strongly nonlinear and, for the two-dimensional case, generalizes 
the classical mean curvature flow equation. We study the motion of one- and two- 
dimensional fronts, finding that the non- homogeneity acts as a "potential function" 
for the motion of the front; i.e., there is wave propagation failure and the steady state 
solution depends on the structure of the function describing the non-homogeneity. 
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1 Introduction 



In this paper we consider the following equation 



<k = D A<p + aP(x,y) [/(</>) + h], (1) 

in a bounded region, Q C R n , n = 1,2, with smooth boundary dQ for Neumann boundary 
conditions on dQ. The positive constants D and a are the diffusion coefficient and the 
production rate of the reactants, respectively. The function / is a bistable function (the 
derivative of a double well potential); i.e., a real odd function with positive maximum equal 
to <f)*, negative minimum equal to —(f)* and precisely three zeros in the closed interval [a_, a+] 
located at a_, ao and a+. For simplicity and without lost of generality we will consider in 
our analysis a_ = —1, ao = and a + = 1. The prototype example is /(</>) = (0 — 3 )/2. 
The constant ft, in (fy), assumed to be small in absolute value, specifies the difference of the 
potential minima of the system as will be explained later. Altough the analysis presented 
below will be valid for a general class of positive different iable function /3, we have in mind 
some particular cases which are described below. In what follows rj is a positive constant. 

Case 1) There is a sequence of points on the real line, k = 1, . . . , N, with iV finite or 
infinite, where the function j3 reaches a maximum, 

N 

(3(x) = ]T e -^- Xk)2 . 

k=l 

Case 2) There is a sequence of lines in the plane, yk, k — 1, . 
where the function j3, independent of x, reaches a maximum, 

JV 

(3{x,y) = Y J e-^y-y^ 2 . (3) 

k=l 

Case 3) There is a sequence of points in the plane, (x^, yj), k — 1, . . . , N, j — 1, . . . , M with 
iV and M finite or infinite, where the function (3 reaches a maximum, 
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(2) 

N, with iV finite or infinite, 



N M 

Pfay) = ££cr(x-x fc , y-yj; rj), where a(x, y; rj) = e~ v{x +y K (4) 

fc=i j=i 

Case 4) There is a sequence of circles in the plane, p — p^, k — 1, . . . , N, with iV finite or 
infinite, and where p represents the radial polar coordinate, where the function (3 reaches a 
maximum, 

N 

=£(*,!/) = £ e~^-^\ (5) 
fc=i 

For sufficiently large values of 77, (3 as given by (@-^D are approximations of distributions 
of discrete sources of reaction. We will refer to the points Xk and (xk,yj),k = 1,...,N, 

j = 1 M as quasi-discrete sources of reaction or quasi-discrete (QD) sites and to the 

stripes y = yu and circles p — p/~ k — 1, . . . , N as quasi-semi-discretes sources of reaction or 
quasi-semi-discrete (QS) sites. We define d to be the shortest distance between two adjacent 
QD or QS sites. 

In the last several years, partial differential equations with nonlinear discrete sources of 
reaction (NDSR) have been used to model phenomena in different fields ranging from physics 
to biology, including the study of pinning in the dislocation motion in crystals, breathers 
in nonlinear crystal lattices, Josephson junction arrays and the biophysical description of 
calcium release waves [?]-[?]. The non-homogeneous version of the bistable equation with 
discrete sources of reaction 

4> t = D A0 + a S(x - x k ) [/(</>) + h], (6) 

k 

has received special attention for its applicability to the dynamics of charge density waves [?] , 
[?]-[?], [?, ?] and to the dynamics of calcium release waves [?, ?, ?]. Equation @ describes the 
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evolution of some concentration (in chemical or biological applications) or order parameter 
(in some physical applications) in a discrete array of nonlinear reaction sites embedded in a 
continuum. In [?] equation (|l|) with an additional term on the right side, —a <ft, and j3 given 
by (0) has been used to model calcium release and uptake in cardiac cells via ryanodine 
receptors. In this model <fi represents the concentration of Ca 2+ and f(<ft) represents the 
calcium-induced calcium release (CICR) activity of the release mechanism. When a = 1, the 
model allows for continuous spatial uptake, whereas when a = 0, it is assumed that release 
and uptake both occur at the QD or QS sites. In [?] the function /(0) was taken to be the 
derivative of a sine-Gordon potential; i.e., f(<f)) = — sin(</>). Note that this last function is 
equal to the derivative of a double well potential, as described above, in a restricted domain 
of definition. 

When the function J2k $(. x ~ x k) in © is replaced by a constant, say 1, then by appropriate 
rescaling we have the bistable equation 

<l> t = bA<f> + f{<f>)+h, (7) 

which describes a phase transition dynamics process, where is a non-conserved order pa- 
rameter. Note that for the particular cases /(</>) = (<p — <P 3 ) / 2 and f(4>) = sin(<f)) equation 
([?]) is the Ginzburg-Landau equation and the overdamped sine-Gordon equation respectively. 
Equation (|7|) can be derived by considering a physical system whose free energy is assumed 
to be of the form 

FM) = In (} (V0)2 + ~ h *) d% ' (8) 
where F is a double well potential having the two equal minima. Note that F(<f) — h <f> is 
a double-well potential with one local minimum and one global minimum. The functional 
derivative of (§) is given by 
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-^ = -bA<t>- m - h, (9) 

where = —F'((f>). The right side of equation (j^j may be considered as a general- 

ized force indicative of the tendency of the free energy to decay towards a minimum. The 
bistable equation (|7|) is obtained by assuming that decreases at a rate proportional to that 
generalized force. 

Equation (|7]) for n = 1 possesses a travelling kink solution moving with velocity proportional 
to h. A kink is a solution that connects the two local minima of the double well potential. 
If h ^ then the kink propagates from the locally stable minimum to the globally stable 
minimum. Of special interest are kinks in which the transition between the two minima takes 
place in a region of order of magnitude e < 1; i.e., the kinks have rapid spatial variation 
between the two ground states. For this case, the point on the line (for n = 1) or the set 
of points in the plane (for n = 2) for which the order parameter vanishes are called the 
interface or the front. Allen and Cahn [?] and Rubinstein, Sternberg and Keller [?] showed 
that for (0) with b <C 1 and h = curved fronts in the plane move with normal velocity 
proportional to their curvature, according to the FMC (flow by mean curvature) equation 



h, (10) 



where y = s(x,t) is the Cartesian description of the interface in the plane, and h is pro- 
portional to h (see also [?]). For a circular interface and h = (both phases have equal 
potential), the curvature is the reciprocal of the radius R, and (|T0| ) becomes R t = — 



whose solution satisfying -R(O) = R is given by R(t) = \j Rq — 2t; i.e., circles shrink to a 

R 2 ~ 

point at a critical time t Cj o = For h > the value of the critical time decreases; i.e., 
tch < tc,o- F° r h < there exists a critical value h c such that if h > h c , then circles still 
shrink to a point in a finite time t c ^ < t Cj o, whereas if h < h c , then circles grow unbound- 
edly. For more general shapes, such an expression for the distance of every point in the 



interface from the origin is difficult to obtain, but there are some analytical results showing 
that the behaviour is similar. Gage and Hamilton [?] proved that (|T0D shrinks convex curves 
embedded in the plane R 2 to a point. They showed that such curves remain convex and 
become asymptotically circular as they shrink. Grayson [?] extended this result to a more 
general case showing that embedded curves become convex without developing singularities; 
i.e., curve shortening shrinks embedded plane curves smoothly to points, with round limiting 
shape. 

On rescaling (|6]) x — > x/d and t —>■ a t and defining b = D/a d 2 (for n = 1 ) [?], equation 
(|) becomes 

(j>t = b A0 + ]T<Kx-x fc ) \f(<f>) + h]. (11) 

k 

The parameter b can be thought of as a measure of how close (|IID is to its continuous limit 
(0). If b —>■ oo equation (|TT|) behaves like its continuous counterpart (with the corresponding 
rescaling); i.e., it possesses a travelling kink solution moving with velocity proportional to h 
[?]. For b small, Mitkov et al. [?, ?] have found numerically that the front dynamics results 
in burst waves characterized by time periodicity in a frame moving along with the front. For 
b small enough, the wave no longer propagates, but relaxes to a stationary kink; i.e., the 
waves are pinned. 

For the one- dimensional version of ([I]), with /3 given by fl2|), as well as for the continuous 
spatial uptake version, a = 1, Keener [?] demonstrated the failure of wavefront propagation 
if the separation between QD sites is large enough. 

For ([[]) we define the spatial and temporal dimensionless variables as in [?] 

x y a 

x = -, V=j, t = at, (12) 

a d 

and we also define the following dimensionless parameters 
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- 1 D 

d V a ' 



rj = rj d 2 



e 



(13) 



Substituting flT2] ) and (T3) into (|T|), dropping the " from the variables and parameters and 
further rescaling the time variable we obtain 



e 2 i = 6 2 A<j> + f3(x,y) [/(0) + e/ i ], 



(14) 



We will consider the case 0< e< 1; i.e., when diffusion is slow, d is large or reaction is fast. 



In Section 2 we make a formal asymptotic analysis to derive an equation of motion for the 
front in equation (0), which in Cartesian coordinates reads 



i - 



___ Sxx S x P x (x, s) P y (x, s) , 

St ~TT4 + 2p( x ,s) -2W^7)-^ S) h (15) 

where the parameter h, which is proportional to h, is defined later. Note that equation ( |To] ) 
expressed in polar coordinates reads 

_ pee p-2 pi- p 2 p e f3 e {p,0) (3 r {p,0) R( fl ,i r rifi . 

This equation generalizes the FMC equation (|l0|) with a strong nonlinearity accounting for 
the influence of the function j3 on the front motion. The method we use is the same as that 
used in [?] for the study of the evolution of kinks in the nonlinear wave equation. We present 
it here in some detail for the sake of completeness. In Section 3 we study the evolution of 
one-dimensional fronts by means of (|T5|). We show that for h = the function (3 acts 
as a "potential function" for the motion of the front; i.e., a front initially placed between 
two maxima of (3 asymptotically approaches the intervening minimum, unlike the classical 
homogeneous equation (0), for which fronts, whose motion is governed by (|10D, move with 
a velocity proportional to h. This result is a consequence of the non-homogeneity of the 
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nonlinear reaction term. In Section 4 we study the evolution of two-dimensional fronts by 
means of (|15|). We show analytically that a radially symmetric and non-constant function 
(3 stabilizes a circular domain of one phase inside the other phase analogous to the one- 
dimensional case. This behavior, also arised as a consequence of the non-homogeneity of 



the nonlinear reaction term. The evolution of closed curves according to (|T^) for (3 given 
by (f|) is studied numerically. We observe that closed convex curves evolve to a final shape 
determined by (3. Our conclusions appear in Section 5. 

2 Asymptotic Analysis: Derivation of the Equation of 
Front Motion 

We assume that for small e > and all t G [0, T], the domain Q can be divided into two open 
regions e) and e) by a curve T(t; e), which does not intersect dQ. This interface, 
defined by 



r(t;e) := {x G Q : <f>(x,t;e) = 0} , (17) 

is assumed to be smooth, which implies that its curvature and its velocity are bounded 
independently of e. We also assume that there exists a solution 4>(x,t; e) of ([[]), defined for 
small e, for all x G Q and for all t G [0, T] with an internal layer. As e — > this solution 
is assumed to vary continuously through the interface, taken the value 1 when x G fl+(t; e), 
— 1 when x G Q-(t, e), and varying rapidly but smoothly through the interface. By carrying 
out a singular perturbation analysis for e < 1, we obtain the law of motion of the interface, 
treating it as a moving internal layer of width 0(e). We focus on the dynamics of the fully 
developed layer, and not on the process by which it was generated. 

In Cartesian coordinates the interface is represented by y = s(x, t, e) for e sufficiently small. 
We assume that the curvature of of the front is small compared to its width and define, in 
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a neighborhood of the interface, a new variable 

y - s(x,t,€) 



z 



e 

which is 0(1) as e — > 0. We call $ the asymptotic form of as e — > with 2 fixed; i.e., 

= $(z,x,t,e). (18) 
The field equation ( |H|) in (2, x, t) coordinates becomes 

e 2 $t - e s t $ z = e 2 $ IX -2e s x <5> zx + (1 + s 2 )$ zz - 

-es« + /?(x,s + e^) /(0) + e/3(x,s + ez) /i. (19) 
The asymptotic expansions of $ and S are assumed to have the form 

$ ~ $o + e $1 + c»( e 2) ; as e _, 0. 

Thus 

/3(x,s + ez) = /?(x, s) + e z (3 y (x,s) + 0(e 2 ). 

Substituting into ( |19"D and equating coefficients of the corresponding powers of e we obtain 
the following problems for 0(1) and 0(e) respectively: 

(l + 4)&„ + P(x,8)f($°) = 0, (20) 



( 1 + si ) $L + P{x, s) f($°) $ 1 = (s xx - 

+ 2 s x C - ^(x, s) z /($°) - /?(x, s) (21) 
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In order to solve (Ell we define a new variable 



(22) 



In terms of £, equation (P0|) reads 



+ /(*°) = 0, (23) 

whose solution is $° = the unique solution of = 0, *(±oo) = ±1, ^(0) = 0. 

Thus 



o I P{ x , s )' 



- z 



In terms of £, x and t, equation (|21| ) reads 



(24) 



$1 £ + /'($°) $ x 



+2 s a 



AeQe, g) + /^(x, s) s x 



/9(z,s)a 

It is straightforward to check that satisfies the homogeneous equation 



(25) 



$^ + "I* 1 = 0. 

That means that the operator A defined as follows 



(26) 



A := i- + /'($°) 



(27) 
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has a simple eigenvalue at the origin with as the corresponding eigenfunction. Then the 
solvability condition for the equation (|25|) gives 



^— ^ r /°V') 2 ^+ 

/3{X, s)2 (1 + sl)2 J-oo 



2 S n 



P X {X,S) + Py{x,S) S x 

[2P{x,s)% (1 + sl) 1 * /3(x, 8 )t (1 + 4)1 J 



/oo 



A simple calculation shows that 



/oo /*oo 
-oo J — oo 



(28) 



/OO 1 /'OO /"CO 1 /"CO 

W = / (*') 2 # and / C d£ = - / (*') 2 de 

-oo Z J —oo J —oo £ J —oo 



(29) 



We define 

_ v]/(+oo)-vfr(-oo) , , 

/_« OT »de ' (30) 

Substituting ( |29| ) and ( |50|) into ( pE|) and rearranging terms we get (P~5[) . Note that for 
/(</>) = (Ginzburg-Landau theory), = tanh | and h = 3 h whereas for /(0) = sin0 
(overdamped sine-Gordon), =4 tan _1 e^ — 7r and h = ^ h. 



3 Front Motion in ID 



For a one-dimensional system, equation ( 1~5[) reads 



St 



2/?(s) 



/3*(s) fe. 



(31) 



We will concentrate on functions (3 of the form 
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N 

^2 



M = E e-^y, (32) 

k=l 

altough the same analysis can be done for a general different iable function. 

In order to analyze the motion of the front we need to look at the roots of the function 

g(s) = -P'(s)-2pl(s)h, (33) 

which are the equilibrium points of the interface. As an example, in Figure |l| we can see the 
graph of /3(s) and g(s) respectively for 77 = 1000 and x\ — 2, x% = 1, x 3 = and x 4 = —1, 
x 5 = —2 and various values of h. For h = 0, g(s) has 9 roots in the range considered. Five 
of them e., they correspond to the maxima of (3(s). The other four 

correspond to the minima of /3(s): s\, S2, S3 and S4 from left to right. It is easy to see that 
Sk, k = 1, ... ,4 are stable whereas Xk, k = 1, ... ,5 are unstable. Thus one dimensional 
fronts initially at non-equilibrium points x move until they reach a stable equilibrium point; 
i.e., a front initially at a point x G {xk,Xk+i) approaches asymptotically Fronts starting 
initially at x > x\ or x < xn will move forever. This behavior is in contrast with the classical 



FMC case (|10|), where one dimensional fronts move only if h ^ 0. In order to understand 
the behaviour of g(s) as h increases above zero we can look at a function /3(s) with a single 
peak at Xi = 0; i.e., (3(s) = e~ vs2 . This function will approximate every (0) if 77 ^> 1, so that 
the influence of peaks on one another is very small. In this case g (s) = 2e~ v s [ris-he*]. 
For h — 0, g(s) vanishes at x = x\ = and it is positive for x > and negative for x < 
(see Fig. |]-b). As /i moves from zero, x, the root of g(s), will be given by the solution 
of r] s — h e _22 2~ = 0, an equation that has a solution as long as r] is sufficiently large and 
h = 0(1). If h > 0, x > 0, g(s) is positive for a; > x and negative for x < x. If h < 
x < 0. We can see the shape of g(s) as ^ increases in Figure |I| In summary, as h increases or 
decreases the behavior of the front is similar to the case h = in contrast with the classical 
FMC case ( ]T0| ) where fronts move with a velocity proportional to the value of h. As an 
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illustration, in Figure 0-a and Q-c we can show the graphs of s t as function of s and of s 
as a function of t, respectively, for 77 = 30 and h = 0,0.5 and 1. In Figure ^-b we can see 
the corresponding graph of (3 as a function of s. We observe that the velocity of the front 
initially increases and then decreases as the front "leaves" the area of the peak of f3(s). As 
h increases, the velocity decreases for a given value of s, as does the asymptotic value of s. 



4 Front Motion in 2D 

In this section we present some analytical and numerical results for the front motion of 
closed curves in the plane according to (|T5j ) . The analysis of front motion in two dimensions 
according fll5"|) with a function f3 of type (||) reduces to the analysis of front motion on a line 
and we shall not consider this case further. 

4.1 Radial Symmetry 

For radially symmetric functions, (3 = (3{p), and initial fronts, equation ([16|) reads 

P 2 (3{p) 

As in the previous Section, the analysis presented here can be performed for a general 
differentiable positive function (3, though here we concentrate on a function (3 of the type 

(!)• 

We begin our analysis for the case N — 1; i.e., (3{p) = e~ T,< " p ~ pl ^' 2 ' , and h = 0. For this case, 
equation becomes 



Pt = — +V(p-Pi)- ( 35 ) 
P 



Equation ( p5| ) has only one equilibrium point, p = yq p\ + \jrfpi + 4 77 J /2 77, which is 
unstable. Note that p — > pi as 77 — > 00. Thus, circles with initial radius po < p shrink 
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to a point in finite time t c whereas circles for which po > p grow unboundedly. This is 
in contrast with the classical case ([R]);i.e., for a constant (5, where circles with any initial 
radius shrink to a point in finite time unless h 7^ (see introduction). For p\ = 0, p(t) = 
y/l/r] + {pi - l/r/) e 2 ^ and then t c = - ln(l - T]pl)/2n. 

For N > 1 we need to look at the roots of the function 

f(p) = -2 (5{p) - p P P {p) - 2 p Pkp) h. (36) 

As an illustration, consider Figure [| which shows the graphs of (3(p) and f(p) for N = 5, 
h = and 77 = 100,30, 11 and 4. In Figures |3|-a and [5|-b, we observe that f(p) vanishes 9 
times; i.e., equation ( |34"D has 9 equilibrium points which are such that f(p) vanishes near 
the maxima or minima of (3{p) with the exception of the first maximum. We call pM,k, 
k = 1, . . . , 5 and p m ,k, k = 1, . . . , 4 the odd and even equilibrium points of (0), respectively, 
starting from that of lower value. We can see that Pm,Hi k = 1, . . . ,5 are unstable whereas 
Pm,k, k = 1,...,4 are stable. Thus circles with an initial radius pM,k < Po < Pu,k+i, 
k — 1, ... ,4 grow or shrink to a circle of radius p m ,k, circles with an initial radius po < Pm,i 
shrink to a point in finite time and circles with an initial radius po > pm,5 grow unboundedly. 
This analysis can be generalized for any value of rj. For a large number of sites pk the analysis 
would be similar. We conclude that the function (3 stabilizes a circular domain of one phase 
inside the other in contrast with the clasical case (|Tt]) where circles of any initial radius 
shrink to a point in finite time. In Figures [5|-c and Q-d we see that for lower values of 77 
some of the equilibrium points disappear (a saddle node bifurcation of stable and unstable 
equilibrium points occurs). For 77 = we expect no equilibrium points, since this corresponds 
to the classical equation (|10"D. In Figure [| we have the graphs of (3{p) and f(p) for N = 5, 
77 = and h = 0, 2 and —2 respectively. We observe the influence of h on the graph of f(p); 
i.e., on the equilibrium points of (|34D and the velocities of circles growing or shrinking. For 
sufficiently large values of h some equilibrium points will eventually dissapear. 



14 



4.2 Some Numerical Results for More General Cases 

In this section we present some numerical results For the evolution of fronts according to 
either ( p~5D or ( |T6D with h = 0. In all cases the function f3(x,y) is given by (f|) with rj = 10, 
N = M = 5 with sites (xk, Vj), k,j = —2, . . . , 2 and (x , yo) = (0, 0). In Figure |5] we show a 
graph of this (3. 

In Figure |6] we see the evolution of a circle of radius 2 for t = 0, 0.15, 0.3. Comparing with 
Figure | we can see that the initial points A = (2, 0), B = (0, 2), C = (-2, 0) and D = (0, -2) 
are relative maxima of (5, while E, F, G and H, the points of intersection between the front 
and the lines y = x and y = —x, lie very near relative minima of f3. As the front evolves, A, 
B, C and D move towards points between the maxima of /3 while E, F, G and H remain 
nearly stationary. The front seeks a position along the minimum of f3. The same behaviour 
can be seen in Figure |7| for the ellipse ^ + y 2 = 1. Because of the very different initial 
conditions, the final front differs front that of Figure ^|. In Figure |8|, we arrive at the same 
final front as in Figure |. In Figure g-a A = (1, 0), B = (0, 2), C = (-1, 0) and D = (0, -2) 
whereas in Figure g-c A = (1.5, 0), B = (0, 1.5), C = (-1.5, 0) and D = (0, -1.5). 



5 Conclusions 



In this manuscript we have derived equation (15) which governs the evolution of a fully 
developed front in a reaction-diffusion system described by (|i~4f) when e< 1 (slow diffusion, 
large separation between sites or fast reaction). This equation generalizes the FMC equation 
(|i0|) to include the effects of stronger nonlinearities and accounts for the influence of the non- 
homogeneous reaction term on the motion of the interface. The motion of fronts according 



to ([151 ) is qualitatively different from that of the homogeneous nonlinear reaction term coun- 
terpart given by the FMC equation, a phenomenon that was pointed out by Keener [?] . This 
difference arises primarily from the fact that the function f3 acts as a "potential" function 
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for the motion of the front. For the one dimensional case, an initial front initially placed 
between two maxima of (3 (which for a homogeneous nonlinear reaction term will move with 
a velocity proportional to h) asymptotically approaches the intervening minimum. For the 
radially symmetric two-dimensional case circular domains of one phase inside the other are 
stabilized. We found numerically that other closed curves present the same phenomenon. 
These results (Fig. |(|||) suggest that the curvature of the front may play a role in "balanc- 
ing" the "force" exerted by the function (3 on the front. Analytical and numerical work will 
be necessary in order to elucidate this effect, which may be crucial in the selection of the 
equilibrium pattern. 



In equation (0) or (|H]) the function (3 can be choosen to depend not only on the spatial 
variable but also on t. The fire-diffuse-fire model of dynamics of intracellular calcium waves 
[?, ?] is of this type. In this case equation fll^D will still govern the evolution of fronts, where 
now (3 = /3(x,y,t). The form of (3(x,y,t) will depend, of course, on the particular model. 
One might, for example, have the product of a spatially dependent function (3(x,y) with a 
probabilistic time dependent function. 

The results presented here have implications for the selection of patterns in systems of the 
type 

' e 2 <j> t = e 2 A<f> + u [f((f)) + eh], 

(37) 

t u t = D u Au + g(u,<f)), 
where D u is a diffusion constant, g is a given nonlinear function the time constant, r, 
is assumed to be large. If u rapidly approaches a non-homogeneous steady state, then 
depending on the initial conditions, it may induce a non-homogeneous steady state in 
where u acts as the function f3 in flip) . 
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Captions 



Figure [I]: 

a) Graph of /3(s) = ELi e- 1000 ^ 2 ^ 2 for x x = 2, x 2 = 1, x 3 = and x A = -1, x 5 = -2. 

b) Graph of g(s) = —(3'(s) — 2 (3^ (s) h for /3(s) as in (a) and h = 0. 

c) Graph of g(s) = — f3'(s) — 2 /33(s) h for as in (a) and h = 10. 

d) Graph of g(s) = -f3'(s) - 2 /3§(s) h for /3(s) as in (a) and h = 20. 
Figure ||: 

a) Dependence of the front velocity s t on the position of the front s in equation ( |5T] ) 

b) Dependence of /3(s) on the position of the front s. 

c) Dependence of the position of the front s on time t. 

In all the cases (3(s) given by fl3|) with N = 2, rj = 30, At = 0.0001 and h = 0, 0.5, 1 from 
above,. The front is initially in so = 0.1. 

Figure |3|: 

Dependence of /3(p) = £Li e"^-^ 2 and /(p) -2 /?(p) -p /? p (p) for Pl = 0, p 2 = 1, p 3 = 2, 
P4 = 3, p 5 = 4 and 

a) 7] = 100. 

b) r? = 30. 

c) r] = 11. 

d) T) = 4. 

Figure f|: 

Dependence of /3(p) = £|=i e^"^) 2 and /(p) - 2 /3(p) - p /^(p) -2 P pl(p)h for Pl = 0, 
p 2 = 1, p 3 = 2, p 4 = 3, p 5 = 4, r? = 30 and 

a) h = 

b) /i = 2 

c) /i = -2 
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Figure |5] 

Graph of (3(x,y) = E|=i T, 5 j=i a ( x - x k ,y - where = e~ 10 ( z2 +^ 2 ) with (x k ,yj), 

k,j = -2, ...,2 and (a? ,2/b) = (0,0). 

Figure || 

Graphs of the evolution of a circle with initial radius equal to 2 according to with 
ft, = and P(x, y) = ELi E|=i c(z ~x k ,y - Vj) where £7 (re, y) = e" 10 ( x2 +^) with (x fc , 
k,j = -2,...,2 and (x ,2/o) = (0,0), for 

a) £ = 

b) £ = 0.15 

c) £ = 0.3 

Figure |7| 

Graphs of the evolution of an ellipse, +y 2 = 1 according to (p~5|) with ft = and /3(rr, y) = 
ELi E|=iO"(^ - £jfe,2/ - %) where cx(x, y) = e~ 10 C^+f 2 ) with (x k ,yj), k,j = -2,..., 2 and 
(£0,2/0) = (0,0), for 

a) £ = 

b) £ = 0.15 

c) £ = 0.4 

Figure || 

Graphs of the evolution of the function (cos(6),2sin(6)) according to fllS]) with ft = 
P(x,v) = ELi E^i^ - Xk,V - Vj) where a(x,y) = e" 10 (a;2+?/2) with (x fc ,%), fc, j = 
-2, . . . , 2 and (x , y ) = (0, 0), for 

a) £ = 

b) £ = 0.15 

c) £ = 0.6 
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